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Abstract 

The role of water ice in the solar system is reviewed from a fluid- dynamical point 
of view. On Earth and Mars, water ice forms ice sheets, ice caps and glaciers at 
the surface, which show glacial flow under their own weight. By contrast, water 
ice is a major constituent of the bulk volume of the icy satellites in the outer solar 
system, and ice flow can occur as thermal convection. The rheology of polycrystalline 
aggregates of ordinary, hexagonal ice Ih is described by a power law, different forms 
of which are discussed. The temperature dependence of the ice viscosity follows an 
Arrhenius law. Therefore, the flow of ice in a planetary environment constitutes a 
thermo-mechanically coupled problem; its model equations are obtained by inserting 
the flow law and the thermodynamic material equations in the balance laws of mass, 
momentum and energy. As an example of gravity-driven flow, the polar caps of 
Mars are discussed. For the north-polar cap, large-scale flow velocities of the order 
of 0.1 ... 1 mm a" 1 are likely, locally enhanced by a factor ten or more in the vicinity 
of surface scarps/troughs. By contrast, the colder south-polar cap is expected to 
be almost stagnant. Tidally heated convection is discussed for the example of the 
icy crust of Europa, where a two-dimensional model predicts the formation of an 
upper, conductive lid and a lower, convective layer with flow velocities of the order 
of 100 mm a . Very little is known about the fluid-dynamical relevance of high- 
pressure phases of water ice as well as ices made up of other materials. 



1 Introduction 

Water ice is an abundant material in the solar system. In the inner solar system, it 
plays an important role on Earth and on Mars, where it forms a cryosphere as an active, 
dynamic part of the respective climate system. Approximately 10% of the land surface of 
the present-day Earth, or 14.6 x 10 6 km 2 , are covered by ice sheets, ice caps and glaciers, 
the total volume of which is approximately 28.7 x 10 6 km 3 [I]. By far the largest single 
ice body is the Antarctic ice sheet, which alone contains about 90% of this ice. All of 
these terrestrial ice bodies are subject to gravity- driven glacial flow with typical velocities 
of tens to hundreds of meters per year. Further components of the terrestrial cryosphere 
are the floating ice shelves and sea ice, seasonal snow and ground ice (permafrost). 
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The ice sheets, ice shelves, ice caps and glaciers have formed by accumulated snow- 
fall over centuries, millenia and more. Over these time-scales, the Earth's climate has 
experienced significant changes known as glacial-interglacial cycles, which are driven by 
periodic changes of the orbital parameters obliquity (axial tilt), eccentricity and precession 
( "Milankovitch cycles"). The last 800,000 years have been characterized by a strong dom- 
inance of the 100,000-year eccentricity cycle [29], with a sequence of shorter interglacials 
(warm periods) like the current Holocene and longer glacials (ice ages). At the last glacial 
maximum approximately 20,000 years ago, large parts of north America and Eurasia were 
covered by ice sheets which no longer exist, and the global ice volume was about three 
times larger than at present. 

On Mars, the polar ice caps are one of the most prominent surface features. The 
seasonal caps, which can extend down to latitudes of approximately 55°N/S, consist of 
only some ten centimeters of CO2 snow which sublimes into the atmosphere during the 
respective spring season. By contrast, the smaller residual caps poleward of approximately 
80°N/S survive the summer seasons, and they are underlain by massive topographic struc- 
tures, which are known as the polar layered deposits [31]. The residual caps and the 
underlying layered deposits are considered to be geomorphological units and shall be re- 
ferred to as the north- and south-polar cap (NPC/SPC), respectively. The Mars Orbiter 
Laser Altimeter (MOLA) measurements of the Mars Global Surveyor (MGS) spacecraft 
have provided a very precise mapping of the surface topographies of the polar caps [32| |4"0] . 
Combined with the estimated cap margins and equilibrated ground topographies discussed 
later (Sect. 4.1), this yields for the NPC a volume of about 1.2 x 10 6 km 3 and an area of 
1.1 x 10 6 km 2 , and for the SPC a volume of 1.8 x 10 6 km 3 and an area of 1.7 x 10 6 km 2 . Due 
to isostatic deflection of the underlying lithosphere, the real volumes may be up to 30% 
larger. In any case, the NPC and SPC are the largest known water reservoirs on Mars. 
Further constituents may be dust, CO2 ice and CO2 clathrate hydrate. Comparable or even 
larger amounts of water may be stored as permafrost in the ground [5]; an idea which was 
corroborated by the spectrometric detection of mid- and high-latitude subsurface layers 
enriched in hydrogen, interpreted as ground ice [2]. 

Similar to the situation on Earth, the Martian polar caps are active components of 
the climate system which interact with the atmosphere thermally, orographically and by 
condensation and sublimation processes of water vapour. Their present topographies are 
the result of the climatic history over at least the last millions of years, which were probably 
characterized by climate cycles as a consequence of strong, quasi-periodic variations of the 
orbital parameters obliquity, eccentricity and precession on time-scales of 10 5 -10 6 years 
[19] . This idea is supported by the light-dark layered deposits of the polar caps, which are 
exposed in the scarps and troughs in the ice surface and close to the margins, and which 
indicate a strongly varying dust content of the ice due to varying atmospheric conditions 
in the past. Further, Head et al. [T5] suggested that Mars underwent "ice ages" during 
periods of high obliquity like that from about 2.1 to 0.4 million years ago (with obliquity 
maxima of ~ 35°). These ice ages are supposedly characterized by warmer polar climates, 
enhanced mass loss of the polar caps due to sublimation and the formation of meters-thick 
ice deposits equatorward to approximately 30°N/S. 

For the Moon and Mercury, which are both devoid of any significant atmosphere and 
therefore subject to intensive solar radiation, water ice may be cold-trapped inside per- 
manently shadowed craters at the poles, which is possible due to the very small tilts of 



2 



the rotational axes. Some evidence for this idea was provided by the space-probe Lunar 
Prospector, which detected hydrogen in the vicinity of the lunar poles by neutron spec- 
trometry [7] , and, for the case of Mercury, by terrestrial radar mapping which revealed a 
highly reflective region on the north pole |31j . 

In the outer solar system (beyond the asteroid belt), owing to the very low tempera- 
tures, evaporation of water ice at the surface of planetary bodies becomes so small that this 
substance survives for time-spans comparable to the age of the solar system [SH]- There- 
fore, while in the inner solar system the geological evolution of the planets is dominated 
by rock material, the history of a large number of bodies in the outer solar system is dom- 
inated by water ice and other moderately volatile substances. Of course, this statement 
does not hold for the giant planets, which do not have a solid surface and are made up of 
highly volatile substances, mainly hydrogen and helium. Thus, it refers to the large moons 
of the Jovian, Saturnian, Uranian and Neptunian systems, Pluto and smaller moons, as- 
teroids and planetary ring particles. Of particular interest are the Jovian moon Europa, 
for which there is strong evidence for an outer icy shell underlain by a deep ocean |18j . 
and the Saturnian moon Titan, for which the findings of the recently landed space-probe 
Huygens suggest a solid surface made up mainly of water ice and a methane hydrosphere 
with rainfall, rivers and seas (see http://huygens.esa.int). 

In this review paper, the focus will be on flowing water ice, which can be driven either 
by gravity forces or by thermal convection. The layout is as follows. In Sect. |2j the 
rheology of polycrystalline aggregates of water ice is discussed, including the effects of 
different creep mechanisms, partial melting and impurities (dust), and suitable flow laws 
(stress-strain-rate relations) are given. This allows the formulation of a closed set of model 
equations for the thermo-mechanically coupled problem of ice flow, which is done in Sect. [3j 
Applications to the Martian polar caps are presented in Sect. |4j New results are shown 
for their large-scale dynamics under assumed steady-state conditions (Sect. 4.1), and the 
local influence of the prominent scarps and troughs in the surface is treated (Sect. 4.2). 
Simulations on tidally heated convection in the icy shell of Europa are discussed in Sect. [5j 
Finally, Sect. [61 concludes the paper. 



2 Ice rheology 
2.1 Structure of ice 

Water ice can exist in a great variety of different phases (Fig. [Tj Ref. [2S])- Most of 
these phases form under high pressure, which results in a denser packing of the water 
molecules compared to the "ordinary" ice Ih. The stability range of the latter is for 
pressures P < 200 MPa, which is equivalent to the hydrostatic pressure of an ice layer 
(density p x = 910 kg m -3 ) of approximately 22 km on Earth (gravity acceleration g = 
9.81ms™ 2 ), 60 km on Mars (g = 3.72 ms -2 ) and 165 km on Europa (g = 1.32 ms~ 2 ) and 
Titan (g = 1.35 ms -2 ). We will limit our discussions to the outer shells of the planetary 
bodies, so that only ice Ih needs to be considered in the following. 

Ice Ih forms hexagonal crystals, that is, the water molecules are arranged in layers of 
hexagonal rings (Fig. |2j Ref. [25]). The plane of such a layer is called the basal plane, 
which actually consists of two planes shifted slightly (by 0.0923 nm) against each other. 
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Figure 1: Phase diagram of the ice- water system. Only stable phases are shown. Figure 
by Petrenko and Whitworth [27?| their Fig. 11.2 on p. 253]. 



The direction perpendicular to the basal planes is the optic axis or c-axis, and the distance 
between two subsequent basal planes is 0.276 nm. This leads to the very low packing factor 
of 34%, which is responsible for the density anomaly of ice Ih (smaller density than liquid 
water). Further, the basal planes can glide on each other when a shear stress is applied, 
comparable to the deformation of a deck of cards. This effect is strongly enhanced by 
the existence and generation of dislocations (structural defects) in real crystals, and the 
mechanism is consequently called dislocation creep. 




0.4523 nm 



Figure 2: Structure of an ice crystal. The circles denote the oxygen atoms of the H 2 
molecules, (a) Projection on the basal plane, (b) Projection on plane indicated by the 
broken line in (a). Figure by Paterson [23| his Fig. 5.1 on p. 80]. 

On the macro-scale, ice aggregates on planetary bodies (grounded ice sheet/cap/glacier, 
ground ice, floating ice shell etc.) are composed of a vast number of individual crystals. 
For instance, for terrestrial ice sheets and glaciers, the typical grain size is of the order of 
millimeters to centimeters. Such a compound is called poly crystalline ice. At the time of 
formation, it can be assumed that the orientation of the crystals is completely at random, 
so that the macroscopic behaviour of the compound will be isotropic. Ice-core studies on 
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Earth have revealed that in the course of time anisotropic fabrics can develop due to the 
strain history which a piece of ice experiences during its motion in an ice sheet or glacier 
(e.g. [TJ[35]; see also [28, this volume]). However, since essentially nothing is known about 
the fabric of extraterrestrial ice bodies, we will assume isotropic conditions for simplicity 
in the following. 



2.2 Pressure melting point 

The melting temperature of pure ice Ih depends on the pressure P. It can therefore be 
written as 

T m = T -f(P), (1) 

where T = 273.16 K is the melting temperature for P = 611.657 Pa, that is, at the triple 
point of water. According to [38], the melting-point depression f(P) is in implicit form 



P 



1 - 626,000 (1 - 0~ 3 ) + 197,135 (1 -0 



21. 2\ 







T m = T - f(P) 
T ~~ T 



(2) 



where P < P < 209.9 MPa and T > T m > 251.165 K (Fig. |3j). For moderate pressures, 
Eq. ([2| can be linearized as 

f(P) = PP, (3) 

where (3 = 7.42 x 10 _8 KPa _1 is the Clausius-Clapeyron constant for pure ice [25J. The 
error of this linearization is negligible for P < 10 MPa, less than 10% for P < 50 MPa and 
approximately 30% for P « 200 MPa (Fig. g). 
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Figure 3: Pressure melting point of ice Ih by Wagner et al. [38J [Eq. (|2j) , solid line] and 
linearized relation ^ (dashed line). 

Impurities of different kinds alter the melting point further. For terrestrial, air-saturated 
glacier ice, Paterson [25] reports a modified Clausius-Clapeyron constant of (3 — 9.8 x 
10 _8 KPa _1 for the linearized melting-point depression If salts are present, an ad- 
ditional melting-point depression occurs which depends on the type of salts and their 
concentration. Intermixtures of ammonia, which may play a role in a supposed subsurface 
ocean on Titan below an outer ice-Ih layer [33j and references therein], have a similar 
effect. 
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2.3 Flow of polycrystalline ice 



When a specimen of polycrystalline ice is subjected to a constant normal or shear stress, it 
responds with a permanent deformation, which continues as long as the stress is applied. 
Typically, an initial, instantaneous elastic deformation of the polycrystalline aggregate 
is followed by a phase called primary creep during which the strain rate decreases con- 
tinuously. This behaviour is related to the increasing geometric incompatibilities of the 
deforming single crystals with different orientations. After some time, a minimum strain 
rate is reached which remains constant in the following, so that the strain increases linearly 
with time. This phase is known as secondary creep. Especially in case of high temperatures 
(> — 10°C), at a later stage dynamic recrystallisation (nucleation and growth of crystals 
which are favourably oriented for deformation) sets in, which leads to accelerated creep 
and finally a constant strain rate significantly larger than that of the secondary creep. 
This is called tertiary creep (Fig. |4j Ref. [25]). Therefore, it can be assumed in good 
approximation that in deforming ice masses secondary creep prevails for low temperatures 
(< — 10°C), whereas tertiary creep prevails for high temperatures (> — 10°C), so that the 
strain rate can be expressed as a unique function of the stress, the ice temperature and 
the pressure. 




secondary creep 

► 



Time 



Figure 4: Creep response of a sample of polycrystalline ice to a constant stress for high 
(> — 10°C, marked by Tf) and low (< — 10°C, marked by T|) temperatures. 

On the basis of these considerations, several forms of a non-linear viscous rheology for 
the flow of polycrystalline ice for different stress, strain-rate and temperature regimes have 
been proposed. They have in common to relate the strain- rate tensor D = symgradi; 
(velocity v) to the Cauchy stress deviator t D , and can be subsumed as 

D = EA(T,P) a ^t D , (4) 

where a = [tr (t D ) 2 /2] 1//2 is the effective stress, n is the stress exponent, d is the grain size 
and p is the grain-size exponent (e.g. (6J [TO], [121 E3 [37]). The flow rate factor A(T,P) 
depends via the Arrhenius law 

A(T,P) = A e-( Q+py V RT (5) 

on the absolute temperature T and the pressure P, where A is the pre-exponential con- 
stant, Q is the activation energy, V is the activation volume and R = 8.314 J mol -1 K _1 is 
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the universal gas constant. The flow enhancement factor E is equal to unity for pure ice 
and can deviate from unity due to the softening or hardening effect of impurities in the ice. 
Since poly crystalline ice is described as a density-preserving (incompressible) medium, the 
pressure is a free field (not governed by a material equation), and the full Cauchy stress 
tensor t is related to the traceless deviator t° by 

t = -Pl + t D , (6) 

where 1 denotes the unity tensor. 

Since appropriate values for the activation volume V are poorly constrained and the 
pressure effect is very small for typical thicknesses of ice sheets and caps, we account for 
it in an approximate way by setting V = and measuring the temperature relative to the 



pressure melting point T m (see Sect. 2.2) instead. To this end, the homologous temperature 

T' = T-T m + T = T + f(P) (7) 

is introduced, and the rate factor can be simplified as 

A(T') = A Q e~ Q ' RT \ (8) 

which depends now exclusively on T' [231 137] . For more details on this approach see 
Appendix |A| 

The general power law Q can be inverted as follows. For the effective strain rate 
5 = {tiD 2 /2) 1 ' 2 we obtain 

5 = EA(T)^, (9) 

or equivalently, solved for a, 

a = [EA{T')]- 1/n d p/n 5 1/n . (10) 
Inserting this in Q and solving for t° yields 

t° = [ J BA(T / )]- 1 #([ J BA(T0]~ 1/n # /ri 5 1/n ) 1_n £> 

, rjP/ n 

= [EAiT'^^z^D. (11) 

By introducing the stress enhancement factor E s = E~ x l n and the associated rate factor 
B(T') = [A{T')\~ l l n , this can be written as 

lP = E B B{r) J[ z&D. (12) 

For terrestrial ice, the well-established Glen's flow law [which actually goes back to [23] 
in the general tensorial form] uses the stress exponent n = 3, the grain-size exponent p = 
and for the temperature range X" < 263 K the pre-exponential constant A = 3.985 x 
1CT 13 s _1 Pa -3 and the activation energy Q = 60 kJ mol -1 [2S]. The rheology defined by 
these parameters describes the grain-size-independent flow mechanism of dislocation creep, 
which prevails in terrestrial glaciers and ice sheets. The flow enhancement factor for ice 
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formed during glacial periods is often set to E = 3, interpreted as the softening influence of 
very small amounts of fine dust, approximately 1 mg kg -1 with particle sizes of 0.1 to 2 /mi 
[H] . This softening is attributed to thin films of liquid water which form around the dust 
particles and lubricate ice deformation. However, at the low temperatures expected in the 
NPC this effect will not be present, and direct hardening will be the dominant influence 
of dust (see Sect. 2.5). 



Durham et al. [B] propose an alternative flow law for grain-size-independent dislocation 
creep, based on laboratory creep tests at a confining pressure of 50 MPa. For the temper- 
ature regime T = 195-240 K, which corresponds approximately to T' = 200-245 K, they 
report the parameters n = 4, p = 0, Aq = 1.259 x 10~ 19 s _1 Pa -4 and Q = 61 kJ moP 1 . 

However, for very low temperature and strain-rate conditions, as they are expected 
in the polar ice caps of Mars, it is not clear whether dislocation creep is still the pre- 
dominant creep mechanism in the polycrystalline ice aggregate. There is evidence that 
other, grain-size-dependent flow mechanisms like grain-boundary sliding become favoured 
instead |10l . These can be described by the parameters n = 1.8, p = 1.4, Aq = 6.20 x 



1Q -14 g-l p a -1.8 m 1.4 and q = 4 g k j mol -l ( gee algo j2Jj)_ 

The relative contributions of the several flow laws (Glen, Durham, Goldsby-Kohlstedt; 
the latter will be abbreviated as "GK" in the following) can be estimated as follows. For 
simple shear in the x-z plane and E = 1, Eq. Q reduces to 



7 = 2A{T') 



dp' 



(13) 



where r — t 

the shear rates resulting from Eq. ( 13 ) for T 1 
100 MPa 



u is the shear stress and 7 = dv x /dz = 2D XZ is the shear rate. Figure [5] shows 

200 K and the stress range from 10 kPa to 
1mm and 10 mm have been assumed, based on an 



For GK, the grain sizes d 
estimate for the Martian north polar cap by Greve and Mahajan [12]. This range is also 
typical for terrestrial ice sheets and glaciers. 
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Figure 5: Shear rate 7 vs. shear stress r for simple shear, computed by Eq. (13) for 
V = 200 K. Flow laws: Glen (n = 3), Durham (n = 4), GK (n = 1.8, p = 1.4, d = 1 mm), 
GK (n = 1.8, p = 1.4, d = 10 mm). 



It becomes evident that the relative contributions of the different flow laws vary strongly 
over this stress range. For low stresses, grain-size dependent flow with a low stress exponent 
(GK) dominates, whereas for higher stresses dislocation creep with a higher stress exponent 
(Glen, Durham) becomes more important. For the shown example with T' = 200 K, 
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the crossover stresses are 166.3 kPa for Glen/GK rf=lmm , 11.3 kPa for Glen/GK d=10mm , 
834.0 kPa for Durham/GK d=lmm and 192.7 kPa for Durham / GK^ =1 n mm ■ 

In order to account for the contributions of grain-size-independent dislocation creep 
and other, grain-size-dependent flow mechanisms simultaneously, a modified version of 
the flow law proposed by Pettit and Waddington [27J may be used. It sums up the two 
contributions via 

/ rr™ 1-1 rr n2 ~ 1 \ 

D = [E^T)— + E 2 A 2 (T')— ] j t D , (14) 

where 

A 1 (T') = (A ) 1 e- Q i/ RT ', A 2 {T') = {A Q ) 2 e~^ RT '. (15) 

In this representation, the index "1" refers to the parameters of either Glen's or Durham's 
flow law, the index "2" to those of the Goldsby-Kohlstedt flow law. However, this approach 
will not be further pursued here. 

Note that, by introducing the viscosity r], the flow law Q with the Arrhenius law in 
the form ^ can be written in compact form as 

1 1 d p 
D = — t D , where ri(T', a, d) = — — r . (16) 



Similarly, the inverse flow law (12) reads 



t D = 2r)(T>, 6, d) D, where 7}(T', 6, d) = 



Of course, the two representations of the viscosity in (16) and (17) are equivalent. 

Evidently, for stress exponents n > 1, the viscosities (16), (17) become infinite in 
the low-stress/deformation limit a, 5 — > 0, and there has been a long-standing debate 
whether this is physically acceptable. Also, some experimental evidence seems to indicate 
a Newtonian (n = 1) rheology for this regime, even though these findings have been 
questioned (see the discussion by |25]). There are different possibilities of introducing a 
finite residual viscosity r] res . The effective stress a and the effective strain rate 5 in the 



expressions (16), (17) can be replaced by max(cr, cr res ) and max(<5, 5 res ), respectively, where 
<T re s and S ies are small positive parameters. Alternatively, a small residual viscosity may 
be directly added, 

1 dP 

^ (T ' a ' rf) = 2lI(T0^ + r/res ' (18) 

or 

, . r IN E S B(T') dPl n , . 

Tj(T',5,d) = ^ ^ sl _ 1/n + Vres- (19) 

The disadvantage of these forms is that they are not equivalent anymore, and an inversion 
according to Eqs. ([9j-( 12) can no longer be executed. 



2.4 Influence of liquid water 

If the temperature of ice reaches the pressure melting point, liquid water may be present 
within the ice matrix, which is stored as lenses at the grain boundaries and in capillary 
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channels at triple grain junctions. Let ip w be the volume fraction of water, then the density 
of the ice-water mixture is 

p = (1 - p w )pi + ip w p w , (20) 

where p\ = 910kgm~ 3 is the density of pure ice and p w = 1000 kg m~ 3 the density of 
pure water. Further, it is clear that the presence of liquid water will reduce the viscosity. 
Tobie et al. [36J propose an exponential softening, which can be expressed by a stress 
enhancement factor E B < 1, 

E s = e~ K ^ } (21) 
or, equivalently, by a flow enhancement factor E > 1, 

E = E~ n = e nK ^, (22) 



where b w = 45. This value was chosen such that 5% melt decreases the viscosity in Eq. (17) 
by an order of magnitude, that is, E s (<p w = 0.05) = 0.1. In terrestrial glaciology, a linear 
relation reported by Paterson [23] has been widely used instead, which reads in terms of 
the enhancement factors 

E = l + b wVw & E B = E-V n =(l + b w <p w )- 1 ' n , (23) 

where b w = ^ = 181.25. This relation has been established based on laboratory measure- 
ments for small water contents below 1% (ip w < 0.01). 

2.5 Influence of dust 

Satellite imagery shows that parts of the Martian polar ice caps appear dark, which in- 
dicates that they consist of ice with some amount of mixed-in dust. Greve and Mahajan 
[12] have laid down that this affects the ice flow in a multiple way, in that direct harden- 
ing can be partly compensated or even overcompensated by the increasing density, which 
increases the driving stresses, and the decreasing heat conductivity, which makes the ice 
at depth warmer and therefore softer. Therefore, the average volume fraction (p d of dust 
is introduced, and the density, p, heat conductivity, n, and specific heat, c, of the ice-dust 
mixture are computed as volume-fraction-weighed averages of the values for pure ice and 
crustal material, 

p = (1 - ipd ) Pi + <p A p C: 

K = (1 - (p d )Ki + Lf d K c , (24) 
PC = (1 - fdjPiCi + VdPcCc, 

with the following parameters: ice density p\ = 910kgm~ 3 , heat conductivity of ice 
«i = 9.828 e -o-0057T[K] Wm -i K -i ) specific h eat Q f ice Ci = (146.3 + 7.253 T[K]) J kg^K" 1 , 

density of crustal material (dust) p c = 2900 kgm -3 , heat conductivity of crustal material 
(dust) k c = 2.5 Wm~ 1 K~ 1 , specific heat of crustal material (dust) c c = 1000 Jkg _1 K _1 . 



The additional density factors in Eq. (24 (3 are necessary because the averaging procedure 
requires volumetric quantities, whereas the specific heat is taken per mass unit. 

Direct hardening is described by a stress enhancement factor E s > 1 based on labora- 
tory measurements of the deformation of ice-dust compounds, 

E s = e 6d¥,d , (25) 
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where = 2 and ip^ < 0.56 [6]. This is equivalent to a flow enhancement factor 

E = E~ n = e ~ nbdlpd . (26) 

Hence, for given stress, temperature and grain-size conditions and a stress exponent n = 3 
a dust content of 10% {tp& = 0.1) leads to an almost twice as hard material (E = 0.55) 
compared to pure ice. 



3 Model equations for ice flow 



The thermo-mechanical problem of ice flow in a planetary environment can be described by 
the balance equations of mass, momentum and energy (e.g. [IE]). For a density-preserving 
(incompressible) medium, which holds for ice in good approximation despite the variability 



expressed by Eqs. (20) and (24)i, the mass balance (continuity equation) reads 

divi; = 0. 



(27) 



For evolving ice bodies, it is convenient to vertically integrate the continuity equation. 
This yields the ice-thickness (H) equation 



dH 

~dt 



-div Q + a s 



<2b, 



(2f 



where Q is the volume flux (vertically integrated horizontal velocity), and a s and a b are 
the mass balances at the surface (positive for supply) and the bottom (positive for loss), 



respectively. The momentum balance yields with the flow law (|17|) the Stokes equation 
— grad P + div r\ (grad v + (grad v) T 



+ pg = 



(29) 



(g: vectorial gravity acceleration), in which the acceleration term pdv/dt has been ne- 
glected due to the very low flow velocities to be expected. From the energy balance, 
Fourier's law of heat conduction 



q = —k(T) gradT 
(q: heat flux, k: heat conductivity) and the caloric equation of state 



(30) 



// = / c(T) dT 

T 



(31) 



(u: specific internal energy, c: specific heat), the temperature-evolution equation 

( dT \ 
pc ( + v ■ grad T j = div (k grad T) + At] 5 2 + r 



(32) 



results. In this relation, the production term 4rj5 2 is the strain heating, and the source 
term r denotes the volumetric heating due to radiation and tidal dissipation. 
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The above equations need to be complemented by dynamic and thermodynamic bound- 
ary conditions at the surface and the bottom of the respective ice body. If we assume that 
the surface is in contact with the atmosphere, then it can be described in good approxi- 
mation as stress-free, that is, 

tn\ s = (33) 

(where n is the outer normal unit vector, and the subscript "s" denotes the surface). The 
surface temperature T s can be prescribed directly as a Dirichlet condition. 

If the bottom is a rigid ice/rock, ice/regolith or ice/sediment interface, no-slip condi- 
tions can be employed, 

v h = (34) 

(the subscript "b" stands for the bottom). As for the temperature field, let us assume that 
the basal heat flux into the ice, q^, is known. This yields the Neumann condition 

= <?b, (35) 

b 

where n is again the outer normal unit vector. 

The situation is different if the bottom is an ice/water interface. In this case, the basal 
stress conditions are governed by the hydrostatic pressure Pb of the water at the interface, 

t- n\ h = -P h n, (36) 

and the bottom temperature equals the pressure melting point, 

T b = T m . (37) 

Provided that the role of impurities is negligible, T m can be obtained from Eq. Q. 



K- 



dT 
dn 



4 Polar caps of Mars 
4.1 Large-scale simulations 

The dynamic and thermodynamic state of the present-day polar caps of Mars will now be 
simulated with the ice-sheet model SICOPOLIS ("Simulation COde for POLythermal Ice 
Sheets"). This model was developed in the mid-1990's for terrestrial applications and has 
later been adapted to the north-polar cap of Mars (see [12], and references therein). It 
solves the ice-flow equations described in Sect. [3] based on the shallow-ice approximation 
[TBI 122], that is, the flow regime is assumed to be simple, bed-parallel shear, the pressure 
is hydrostatic and lateral shear stresses as well as normal stress deviators are neglected 
(see also [3j this volume]). Inputs from the environment are specified by the mean annual 
surface temperature, the net surface mass balance (ice accumulation minus ablation) and 
the basal ( "geothermal" ) heat flux from the underlying lithosphere. The numerical solution 
of the model equations is carried out by a finite-difference integration technique. 

Here, we only consider steady-state conditions, that is, the ice surface is held fixed 
as given by the Mars Orbiter Laser Altimeter (MOLA) data of the Mars Global Surveyor 
(MGS) spacecraft [32ll4*0] . This makes it unnecessary to prescribe the surface mass balance. 
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Figure 6: North-polar cap of Mars: MOLA surface topography h |32j, computed ice thick- 
ness H, computed surface velocity v B and computed homologous basal temperature T£ for 
present-day steady-state conditions. Heavy-dashed lines indicate the ice-cap margin. 



The topographies are shown in Figs. [6] and [7] (top left panels) for the north- and south- 
polar caps and their surroundings. It is striking that the overall shape of the north-polar 
cap is quite regular and smooth, whereas the south-polar cap appears much more rugged. 

The ice temperature, velocity and bottom topography are allowed to evolve freely until 
the steady state is reached. Surface temperature is delivered by the model by Grieger 
(pers. comm. 2004; see also |http:/ /www.space-vision.biz/marstemperatures.html[ ), which 
is based on a zonal and daily mean energy balance including a simple scheme for CO2 
condensation and evaporation. Results agree well with those of the Martian Climate 
Database [21]; however, both methods fail in reproducing the observed year-round CO2 
cover of the southern residual cap. Therefore, the southern mean-annual surface temper- 
ature is corrected such that it equals the C0 2 sublimation temperature of — 128°C within 
85°S. Further, the geothermal heat flux is assumed to be 35mWm -2 [30], and the source 



term r in the temperature evolution equation (32) is neglected. 



The unknown topography of the solid ground below the polar caps is computed in 



13 




x (km) x (km) 



Figure 7: South-polar cap of Mars: MOLA surface topography h [32J, computed ice thick- 
ness H, computed surface velocity v B and computed homologous basal temperature T£ for 
present-day steady-state conditions. Heavy-dashed lines indicate the ice-cap margin. 

two steps. First, the equilibrated ground for ice- free conditions, which is required as 
a reference topography, is determined by a smooth extrapolation of the ice-free ground 
surrounding the north and south polar cap, respectively [13]. Second, the actual ground 
topography is obtained by superposing the isostatic deflection (downward displacement w) 
of the underlying lithosphere due to the ice load. By modelling the lithosphere as a thin 
elastic plate, the isostatic deflection is governed by the bi-potential equation 

K x V 4 w = pgH - p a gw, (38) 

where K\ = 10 25 N m is the flexural stiffness of the lithosphere, V 4 the bi-potential operator 
in the horizontal plane and p a = 3500 kg m~ 3 the density of the asthenosphere (viscous 
mantle layer below the elastic lithosphere). The two load terms on the right-hand side 
are the ice load itself {pgH) and the counteracting buoyancy force which the deflected 
lithosphere experiences from the asthenosphere below (p^gw). Note that the ice thickness 
H depends on w and is therefore part of the solution. In transient scenarios, which are 
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not considered here, the elastic deflection (38) is not assumed instantaneously, so that an 
additional evolution equation is required for the non-equilibrium displacement [TTI |2"0] . 
The simulations have been carried out by applying Glen's flow law as discussed in 



Sect. 2.3 The dust content is assumed to be 20% (y?d = 0.2, see Sect. 2.5). Horizontal 
resolution is 20 km in the stereographic plane with standard parallel 71°N/S, vertical 
resolution is 51 grid points in the cold-ice column and 11 grid points in the lithosphere 
column, and the time-step is 1000 years. 

Results for the ice thickness, the surface velocity and the homologous basal temperature 
are shown in Figs. [6] and [7] (top right and bottom panels) for both polar caps. The north- 
polar cap assumes its maximum thickness of 3.62 km almost exactly at the pole, whereas 
the maximum thickness of the south-polar cap, 3.88 km, is offset by approx. 150 km in 
10°W direction. Of the entire ice volume of the NPC, which is 1.53 x 10 6 km 3 , a fraction 
of 1.23 x 10 6 km 3 (80%) is above the level of the equilibrated ground, and the remaining 
20% are due to the isostatic deflection of the lithosphere. This ratio is almost the same 
for the SPC, with a total volume of 2.30 x 10 6 km 3 and a volume above the equilibrated 
ground of 1.85 x 10 6 km 3 . 

Compared to terrestrial ice sheets and glaciers with typical flow velocities of tens to 
hundreds of meters per year, the flow of the Martian polar caps is very slow. The NPC 
reaches a maximum surface velocity of 0.98 mm a -1 , and south of 85°N surface velocities 
are everywhere less than 0.01 mm a" 1 . Therefore, the active, dynamic zone only consists 
of the interior, thick part of the ice cap. The situation is even more extreme for the SPC, 
where the maximum surface velocity is as small as 0.055 mm a" 1 , and by far the largest 
part of the ice cap flows at speeds slower than a micrometer per year (note the different 
scales of the colour bars). Therefore, the SPC can be considered as essentially stagnant. 

Note that the 20-km resolution applied here does not resolve small-scale structures like 
the scarp/trough systems, which may lead to locally enhanced flow velocities. This will be 
discussed below (Sect. 4.2) for the NPC. 

The slow flow velocities of the Martian polar caps lead to a conduction-dominated heat 
transport within the ice. The ratio of convection to conduction is given by the Peclet 
number 

A=™, (39) 
[a] 

where [U] is a velocity scale, [L\ is a scale of extent and [a] is a scale for the thermal 
diffusivity a = K/(pc). With [U] = 0.1 mm a -1 for the NPC and 0.01 mm a -1 for the 
SPC, respectively, [L] = 300 km (only the inner, dynamic region considered; see velocity 
panels of Figs. @ and Q and [a] = 2 x 10- 6 m 2 s -1 (for -100°C and 20% dust), one 
obtains Peclet numbers of Pe ~ 0.5 for the NPC and Pe ~ 0.05 for the SPC, respectively. 
Evidently, both values are less than unity, so that conduction outweighs convection. This 
leads to temperature profiles which increase essentially linearly with depth (not shown). 
By contrast, for terrestrial ice sheets and ice caps, Pe ^> 1 always holds, so that their 
temperature fields are mainly governed by heat convection. 

The different flow behaviour of the NPC and SPC is mainly due to the lower ice 
temperatures of the SPC, which are a consequence of the ~ 20°C difference in surface 
temperatures. Since the thicknesses of both caps are similar, and the temperature dis- 
tribution is mainly controlled by heat conduction, this transfers directly to a ~ 20°C 
difference in average basal temperatures (again, note the different scales of the colour bars 
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in the respective panels of Figs. [6] and [7]). The highest homologous basal temperatures 
are — 69.1°C for the NPC and — 89.7°C for the SPC, respectively, which demonstrates the 
insulating effect of the ice caps against the much colder surface temperatures. Neverthe- 
less, the basal temperatures are in any case far below the pressure melting point. This is 
a very robust result, so that the presence of large amounts of subglacial liquid water as a 
potential habitat for Martian lifeforms can essentially be ruled out. 

The simulation for the north-polar cap has been re-run with the alternative flow laws 
shown in Fig. [5] (GK with d = 1mm and 10mm, Durham). As it was already found 
by Greve and Mahajan [12] . this has a very significant influence on the computed flow 
velocities. The maximum surface velocity varies by almost three orders of magnitude, 
from the largest value 3.02 mm a -1 for GK^ =lmm via the above-mentioned 0.98 mm a" 1 for 
the reference simulation (Glen) and 0.12 mm a -1 for GKd =10mm to only 0.013 mm a -1 for 
Durham. As expected, this order corresponds to the low-stress regime in Fig. [5] Since 
the creep mechanism is probably a combination of dislocation creep (Glen, Durham) and 
grain-size-dependent creep (GK), we consider a value of the order of 1mm a -1 as most 
likely. By contrast, the basal temperature is virtually unaffected by the assumed flow 
law, which is a consequence of the dominance of flow-independent heat conduction over 
flow-dependent heat convection. 



4.2 Detailed simulations of the scarps and troughs 



In the large-scale simulations discussed above (Sect. 4.1), the spiralling scarps and troughs 
which cut up to several hundred meters into the surface of the Martian polar caps, have 
not been resolved. This was tackled in a study by Hvidberg |17J, where a flowline of the 
north-polar cap extending from the pole in 160°E direction was investigated in detail, and 
an axisymmetric ice cap was assumed. The flowline was discretized by a finite-element 
grid at kilometer-scale horizontal resolution, and the ice- flow equations of Sect. [3] were 
solved without further approximations. Otherwise, the model set-up is very similar to 
that described in Sect. 14. 1L 

The temperature and flow fields computed by Hvidberg [TTJ are shown in Fig. [8j While 
the overall values agree very well with those reported above for the large-scale simulations 



(Sect. 4.1| Fig. [6]), it becomes evident that the topographic disturbances imposed by the 



scarps/troughs result in significant disturbances of the ice flow which propagate all the 
way down to the bottom. The local flow accelerations are highlighted in Fig. [9j which 
shows the surface velocities and the surface mass balance required for maintaining the 
steady state. For the most pronounced trough about 220 km away from the pole, the 
horizontal surface velocity is larger than 15 mm a -1 , and even the vertical velocity reaches 
values of almost ±10 mm a -1 , distributed such that the trough would close in the absence 
of any mass exchange. In order to keep the troughs open, this must be balanced by a mass 
exchange at the surface such that ice accumulates outside the troughs, but is removed 
from within them at a rate of some millimeters per year (Fig. [9j bottom panel). This 
exchange pattern was already proposed by Fisher [8j |9] and termed "accublation model" . 
Physical processes behind the "accublation" exchange are likely differential ablation due to 
the albedo contrast (white material outside vs. darker material inside the troughs) and/or 
enhanced wind erosion due to the formation of local turbulences in the troughs. 
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Figure 8: Computed temperature and velocity field of the north polar cap of Mars, for a 
flowline extending southward from the pole in 160°E direction. Figure by Hvidberg [T7l 
her Fig. 2]. 
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Figure 9: MOLA surface topography, computed surface velocity and computed steady- 
state mass balance corresponding to the simulation shown in Fig. |8j Figure by Hvidberg 
PH her Fig. 3]. 
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5 Icy shell of Europa 



While the polar caps of Mars are localized ice masses on the planet's surface which rest 
on solid land, the icy shell of Europa envelops the entire planetary body and is underlain 
most likely by a deep ocean. Also, the supposed thickness within the range of 10-50 km is 
much larger than that of the Martian polar caps. This allows the formation of convection 
cells, so that vertical motion may play a much larger role. 




Figure 10: Model of the radial structure of Europa with a 20-km thick outer ice shell 
underlain by a deep ocean. The ice shell consists of an upper conductive lid and a lower 
convective layer. A possible temperature profile is shown in the inset plot. Figure by Tobie 
et al. [Ml their Fig. 1]. 



Tobie et al. [35] investigated this problem by assuming a 20-km thick ice shell with a 
surface temperature of 100 K and a bottom temperature of 270 K, which is approximately 
equal to the pressure melting point (Fig. 10). The ice-flow equations of Sect. [3] are solved 
by finite differences in a two-dimensional cross-section of 40 km width, and the domain is 
treated as Cartesian. For the viscosity law, a simplified Newtonian viscosity is assumed, 
that is, n — 1 and p = [see Eqs. (16) and (17)], and the dependence on the homologous 
temperature T' is modelled by an activation energy of Q = 50 kJ mol -1 , which is close to 
the value of Goldsby and Kohlstedt [10]. The bottom viscosity is set to 1.5 x 10 14 Pas in 
the reference simulation. The tidal dissipation r in the energy balance (32) is calculated 
by assuming a viscoelastic response of Europa to the tidal forcing, and the parameters of 



the Maxwellian rheology are prescribed variably in the four layers indicated in Fig. 10 



This can lead to partial melting in the lower parts of the ice shell, the effect of which is 
accounted for by Eqs. (20)-(22). In the vertical component of the Stokes equation (29), 
the buoyancy force due the varying ice density, which results from partial melting and 
thermal expansion, is added. 



The resulting fields of ice flow, temperature and tidal heating are shown in Fig. 11 



As already sketched in Fig. 10 , the huge temperature (and therefore viscosity) differences 
lead to the formation of an upper, conductive sublayer of approximately 8 km thickness, 
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Figure 11: Computed flow field (arrows in both panels), temperature (panel a) and tidal 
heating (panel b) in Europa's ice shell for the reference simulation by Tobie et al. [36, their 
Fig. 4]. The solid contour in panel (a) delimits the area where the temperature is at the 
pressure melting point. 



which is essentially rigid, and in which the temperature increases linearly with depth. By 
contrast, in the underlying convective sublayer significant ice flow occurs in the form of 
convection cells, with a maximum horizontal velocity of 0.29 ma" 1 , a maximum upward 
velocity of 0.265 ma -1 and a maximum downward velocity of 0.352 ma" 1 . Therefore, 
vertical temperature gradients are small in the convective sublayer. In the rising plumes, 
the ice temperature reaches the pressure melting point, and consequently partial melting 
occurs. The tidal heating is negligible in the upper, conductive sublayer, but reaches 
significant values up to 2.7 x 10" 6 Wm -3 in the lower parts of the ice shell. It is also 
interesting to note that the convection cells entail large lateral gradients of the temperature 
and the tidal heating. 

Fig. 12 shows the corresponding heat flux at the surface and the bottom of the ice shell 
(solid line in the plots). As it can be expected, the lateral variability of the heat flux at the 
surface is rather small, and the average value is approximately 40mWm -2 . By contrast, 
at the bottom, the convection cells produce a heat-flux pattern highly variable in space 
and time (the latter is not shown in the snapshot figure) with an average of about one 
fourth of the surface heat flux. 

In addition to the reference simulation discussed above, Tobie et al. [36J studied also 



the effects of different bottom viscosities, tidal dissipation rates (dashed line in Fig. 12 



19 



50 
45 
40 

35 
30 



Surface heat flux (mW/m 2 ) 



10 20 30 40 

X (km) 



Bottom heat flux (mW/m 2 ) 




Figure 12: Surface and bottom heat flux corresponding to the simulation shown in Fig. [TT 
(solid line). The dashed line is for the same simulation without tidal heating |36[ their 
Fig. 6]. 



etc.) and ice-shell thicknesses. This shall not be further reported here. 



6 Concluding remarks 

We have seen that water ice plays an important role in the climate systems of Earth and 
Mars. Further, it is a major constituent of the bulk volumes of the satellites of Jupiter, 
Saturn, Uranus and Neptun as well as the planet Pluto. The ice sheets, ice caps and 
glaciers on Earth and Mars show glacial flow, driven by their own weight. By contrast, ice 
flow driven by thermal convection can occur in the crust and interior of the icy satellites, 
depending on the temperature gradient between the warm interior and the cold surface. 
Here, we have limited the discussions to the "ordinary" ice Ih. However, at greater depths 
the occurrence of high-pressure phases can be expected, the rheology of which is poorly 
known. 

Further "ices" made up of other moderately volatile substances are known or supposed 
to occur, for example, on Mars (C0 2 ), Io (S0 2 ), Europa, Ganymede and Callisto (NH 3 ), 
Titan (CH 4 and other carbohydrates), Triton and Pluto (N 2 , CH 4 ). However, from a fluid- 
dynamical point of view, very little is known about their relevance, and it has not been 
possible to assess so far whether these materials are significant parts of flowing ice masses 
somewhere in the solar system. 
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A Pressure dependence of ice flow 

The flow rate factor (|8l) can be written as 



AiT') 




Aoe * P {-RT(l + f3P/T))> (40) 



where the linearized melting-point depression ^ has been used for simplicity. For the 
stability range of ice Ih, P < 200 MPa, we have (3P < 20 K, and therefore pP/T < 1. 
This allows the Taylor approximation 

A(D = A^ P (-U - = A e,p (- UQ - ^)) • (41) 



RT v T 



By comparing this with the flow rate factor in the general form ([5]), the activation volume 
which corresponds to the simplified form (pi) is 



QI3P Q3 , N 

FV = -^- =► V = -^f. (42) 

For Q = 60 kJ mol" 1 and T = 200 K, this yields V = —2.94 x 10~ 5 m 3 mol -1 , which agrees 
well with the (quite uncertain) values reported by Paterson [25], V = — 1.7 x 10 -5 m 3 mol" 1 , 
and by Durham et al. [6], V = —1.3 x 10~ 5 m 3 mol _1 . Therefore, the pressure dependence 
of the flow rate factor is described reasonably well by the simplified form ^ which depends 
on the homologous temperature only. 
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